/*
 * Class:        NormalInverseGaussianProcess
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
package umontreal.ssj.stochprocess;
import umontreal.ssj.rng.*;
import umontreal.ssj.probdist.*;
import umontreal.ssj.randvar.*;
import umontreal.ssj.hups.*;

/**
 * This class represents a normal inverse gaussian process (NIG). It obeys
 * the stochastic differential equation @cite pBAR98a&thinsp;
 * @anchor REF_stochprocess_NormalInverseGaussianProcess_eq_nig
 * @f[
 *   dX(t) = \mu dt + dB(h(t)), \tag{nig}
 * @f]
 * where @f$\{B(t),  t\ge0\}@f$ is a  @ref BrownianMotion with drift
 * @f$\beta@f$ and variance 1, and @f$h(t)@f$ is an
 * @ref InverseGaussianProcess @f$IG(\nu/\gamma,\nu^2)@f$, with @f$\nu=
 * \delta dt@f$ and @f$\gamma= \sqrt{\alpha^2 - \beta^2}@f$.
 *
 * In this class, the process is generated using the sequential technique:
 * @f$X(0)=x_0@f$ and
 * @f[
 *   X(t_j) - X(t_{j-1}) =\mu dt + \beta Y_j + \sqrt{Y_j} Z_j,
 * @f]
 * where @f$Z_j \sim N(0,1)@f$, and @f$Y_j \sim IG(\nu/\gamma,\nu^2)@f$
 * with @f$\nu= \delta(t_j - t_{j-1})@f$.
 *
 * There is one  @ref umontreal.ssj.rng.RandomStream used to generate the
 * @f$Z_j@f$’s and there are one or two streams used to generate the
 * underlying  @ref InverseGaussianProcess, depending on which IG subclass is
 * used.
 *
 * In finance, a NIG process usually means that the log-return is given by a
 * NIG process;  @ref GeometricNormalInverseGaussianProcess should be used in
 * that case.
 *
 * <div class="SSJ-bigskip"></div><div class="SSJ-bigskip"></div>
 */
public class NormalInverseGaussianProcess extends StochasticProcess {
    protected RandomStream streamIG1;   // U[0,1) gen used in the inverse gaussian
    protected RandomStream streamIG2;   // U[0,1) gen used (maybe) in the inverse gaussian
    protected RandomStream streamBrownian;   // U[0,1) gen for the normal "Brownian"

    protected InverseGaussianProcess igProcess; 
    protected NormalGen normalGen;
    // Randomized time generated by the InverseGaussianProcess.
    protected double[] stochTime;
    double[] dt;
    double[] mudt;

    // Parameters of the normal inverse gaussian
    protected double mu;
    protected double delta;
    protected double alpha;
    protected double beta;
    protected double gamma;

   /**
    * Given an  @ref InverseGaussianProcess `igP`, constructs a new
    * `NormalInverseGaussianProcess`. The parameters and observation times
    * of the IG process will be overriden by the parameters of the NIG
    * process. If there are two  @ref umontreal.ssj.rng.RandomStream ’s in
    * the  @ref InverseGaussianProcess, this constructor assumes that both
    * streams have been set to the same stream.
    */
   public NormalInverseGaussianProcess (double x0, double alpha,
                                        double beta, double mu,
                                        double delta,
                                        RandomStream streamBrownian,
                                        InverseGaussianProcess igP) {
        setParams (x0, alpha, beta, mu, delta);
        this.streamBrownian = streamBrownian;
        normalGen = new NormalGen(streamBrownian); // N(0,1)
        igProcess = igP;  // its params will be set in init().
        this.streamIG1 = igProcess.getStream();
        this.streamIG2 = streamIG1;
    }

   /**
    * Constructs a new <tt>NormalInverseGaussianProcess</tt>. The string
    * argument corresponds to the type of underlying
    * @ref InverseGaussianProcess. The choices are SEQUENTIAL_SLOW,
    * SEQUENTIAL_MSH, BRIDGE and PCA, which correspond respectively to
    * @ref InverseGaussianProcess,  @ref InverseGaussianProcessMSH,
    * @ref InverseGaussianProcessBridge and
    * @ref InverseGaussianProcessPCA. The third
    * @ref umontreal.ssj.rng.RandomStream, `streamIG2`, will not be used
    * at all if the SEQUENTIAL_SLOW or PCA methods are chosen.
    */
   public NormalInverseGaussianProcess (double x0, double alpha,
                                        double beta, double mu,
                                        double delta,
                                        RandomStream streamBrownian,
                                        RandomStream streamIG1,
                                        RandomStream streamIG2,
                                        String igType) {
        setParams (x0, alpha, beta, mu, delta);
        this.streamIG1 = streamIG1;
        this.streamIG2 = streamIG2;
        this.streamBrownian = streamBrownian;
        normalGen = new NormalGen(streamBrownian); // N(0,1)

        // The initial time of igProcess here, 0., is arbitrary: 
        // init() sets igProcess.x0 = t[0].
        if(igType.compareTo("SEQUENTIAL_SLOW") == 0)
            // the initial value, 0.0, of the IG will be overriden by
            // the initial time, when it will be set.
            igProcess = new InverseGaussianProcess(0.0, delta, gamma, 
                                streamIG1);
        else if(igType.compareTo("SEQUENTIAL_MSH") == 0)
            igProcess = new InverseGaussianProcessMSH (0.0, delta, gamma, 
                                streamIG1, streamIG2);
        else if(igType.compareTo("BRIDGE") == 0)
            igProcess = new InverseGaussianProcessBridge (0.0, delta, gamma, 
                                streamIG1, streamIG2);
        else if(igType.compareTo("PCA") == 0)
            igProcess = new InverseGaussianProcessPCA (0.0, delta, gamma, 
                                streamIG1);
        else throw new IllegalArgumentException("Unrecognized igType");
    }

   /**
    * Same as above, but all  @ref umontreal.ssj.rng.RandomStream ’s are
    * set to the same stream, `streamAll`.
    */
   public NormalInverseGaussianProcess (double x0, double alpha,
                                        double beta, double mu,
                                        double delta,
                                        RandomStream streamAll,
                                        String igType) {
        this(x0, alpha, beta, mu, delta, streamAll, streamAll, streamAll, igType);
    }

   /**
    * Generates the path. This method samples each stream alternatively,
    * which is useful for quasi-Monte Carlo, where all streams are in fact
    * the same iterator on a  @ref umontreal.ssj.hups.PointSet.
    */
   public double[] generatePath() {
        if(igProcess.getNumberOfRandomStreams() != 1)
            return generatePathTwoIGStreams();

        double x = x0;
        double[] randomNormal = new double[d];
        double[] randomIG1 = new double[d];
        for (int j = 0; j < d; j++)
        {
            randomIG1[j]    = streamIG1.nextDouble();
            randomNormal[j] = streamBrownian.nextDouble();
        }

        stochTime = igProcess.generatePath(randomIG1);
        for (int j = 0; j < d; j++)
        {
            double dy = stochTime[j + 1] - stochTime[j];
            x += mudt[j] + beta*dy + 
                 Math.sqrt(dy) * NormalDistQuick.inverseF01(randomNormal[j]);
            path[j + 1] = x;
        }
        observationIndex = d;
        return path;
    }


    protected double[] generatePathTwoIGStreams() {
        double x = x0;
        double[] uniformNormal = new double[d];
        double[] uniformIG1 = new double[d];
        double[] uniformIG2 = new double[d];

        for (int j = 0; j < d; j++)
        {
            uniformIG1[j]    = streamIG1.nextDouble();
            uniformNormal[j] = streamBrownian.nextDouble();
            uniformIG2[j]    = streamIG2.nextDouble();
        }

        stochTime = igProcess.generatePath(uniformIG1, uniformIG2);
        for (int j = 0; j < d; j++)
        {
            double dy = stochTime[j + 1] - stochTime[j];
            x += mudt[j] + beta*dy + 
                 Math.sqrt(dy) * NormalDistQuick.inverseF01(uniformNormal[j]);
            path[j + 1] = x;
        }
        observationIndex = d;
        return path;
    }

/**
 * Returns the value of the process for the next time step. If the underlying
 * @ref InverseGaussianProcess is of type  @ref InverseGaussianProcessPCA,
 * this method cannot be used. It will work with
 * @ref InverseGaussianProcessBridge, but the return order of the
 * observations is the bridge order.
 */
public double nextObservation() {
       double igNext = igProcess.nextObservation();
       observationIndex = igProcess.getCurrentObservationIndex();
       stochTime[observationIndex]  = igNext;
       double dY = igNext - stochTime[0];
       double dT = t[observationIndex] - t[0];
       path[observationIndex] = x0  +  mu * dT  +  beta * dY  +
                                Math.sqrt(dY) * normalGen.nextDouble();
       return path[observationIndex];
    }


    protected void init()
    {
        super.init();
        igProcess.setParams(delta, gamma);
        if(observationTimesSet){
            stochTime = new double[d+1];
            stochTime[0] = t[0];
            dt = new double[d];
            mudt = new double[d];
            for(int i = 0; i < d; i++){
               dt[i]   = t[i+1] - t[i];
               mudt[i] = dt[i] * mu;
            }
        }
    }

/**
 * Sets the observation times on the NIG process as usual, but also sets the
 * observation times of the underlying  @ref InverseGaussianProcess. It
 * furthermore sets the starting *value* of the  @ref InverseGaussianProcess
 * to `t[0]`.
 */
public void setObservationTimes(double t[], int d) {
        super.setObservationTimes(t, d);
        igProcess.setObservationTimes(t, d);
        igProcess.x0 = t[0];
    }

   /**
    * Sets the parameters. Also, computes @f$\gamma=
    * \sqrt{\alpha^2-\beta^2}@f$.
    */
   public void setParams (double x0, double alpha, double beta,
                          double mu, double delta) {
        // Initializes the parameters of the process
        if (delta <= 0.0)
            throw new IllegalArgumentException ("delta <= 0");
        if (alpha <= 0.0)
            throw new IllegalArgumentException ("alpha <= 0");
        if (Math.abs(beta) >= alpha)
            throw new IllegalArgumentException ("|beta| >= alpha");

        this.x0    = x0;
        this.mu    = mu;
        this.delta = delta;
        this.beta  = beta;
        this.alpha = alpha;

        gamma = Math.sqrt(alpha*alpha - beta*beta);
        if (observationTimesSet) init();
    }

   /**
    * Returns alpha.
    */
   public double getAlpha() {
      return alpha;
   }

   /**
    * Returns beta.
    */
   public double getBeta() {
      return beta;
   }

   /**
    * Returns mu.
    */
   public double getMu() {
      return mu;
   }

   /**
    * Returns delta.
    */
   public double getDelta() {
      return delta;
   }

   /**
    * Returns gamma.
    */
   public double getGamma() {
      return gamma;
   }

   /**
    * Returns the analytic average, which is @f$\mu t + \delta t \beta/
    * \gamma@f$.
    */
   public double getAnalyticAverage (double time) {
        return mu*time+delta*time*beta/gamma;
    }

   /**
    * Returns the analytic variance, which is @f$\delta t \alpha^2 /
    * \gamma^3@f$.
    */
   public double getAnalyticVariance (double time) {
        return delta*time*alpha*alpha/gamma/gamma/gamma;
    }

   /**
    * Only returns the stream if all streams are equal, including the
    * stream(s) in the underlying  @ref InverseGaussianProcess.
    */
   public RandomStream getStream() {
        if( (streamIG1 != streamIG2)                ||
            (streamIG1 != streamBrownian)         ||
            (streamIG1 != normalGen.getStream())  ||
            (streamIG1 != igProcess.getStream())  )
            throw new UnsupportedOperationException
            ("Two different streams or more are present");
        return streamIG1;
    }

   /**
    * Sets all internal streams to `stream`, including the stream(s) of
    * the underlying  @ref InverseGaussianProcess.
    */
   public void setStream (RandomStream stream) {
        streamIG1 = streamIG2 = streamBrownian = stream;
        normalGen.setStream(stream);
        igProcess.setStream(stream);
    }

}